#!/usr/bin/env python3
# -*- coding:utf-8 -*-

'''
author: maoxuner
e-mail: maoxuner@126.com
date: 2016年4月24日
filename: stability_analysis_vb.py
version: 1.1
details: 本程序从VB版本转化过来，详细功能参见VB版本的代码
warning: 您可以任意复制、传播本程序，但是请保留此注释内容。
         如果您要基于此程序进行二次开发，请补充注释修改的内容，
         并把新的程序无偿的分享给大家。谢谢！
'''

from math import sin, cos, tan, atan, sqrt

def get_range(start, end, step):
    tmp = start
    while tmp <= end:
        yield tmp
        tmp += step

ed = 1.0
kc = kc_dz = 5.0

print('欢迎使用本程序计算坝坡稳定，请按照文字提示输入数据。')
h = float(input('请输入坝高：'))
d = float(input('请输入坝宽：'))
m = float(input('请输入平均坡度：'))
wh = float(input('请输入水深：'))
r = float(input('请输入砂砾料天然容重：'))
rf = float(input('请输入砂砾料湿容重：'))
dz = input('请输入地震设计烈度（7~9）：')

factors_quake = {'7': (3.0, 0.1), '8': (2.5, 0.2), '9': (2.0, 0.4)}
am, ah = factors_quake[dz]

xb = h * m
for xa in get_range(0.0, xb - ed, ed):
    ya = xa / m
    for xd in get_range(xa + ed, xb, ed):
        for yd in get_range(ya + ed, xd / m - ed, ed):
            a2 = atan((yd - ya) / (xd - xa))
            ye = xd / m
            xc_max = xd + d
            if xa + (h - ya) / tan(a2) < xc_max:
                xc_max = xa + (h - ya) / tan(a2)
            for xc in get_range(xb, xc_max, ed):
                if xd != xc:
                    a1 = atan((h - yd) / (xc - xd))
                    if wh <= ya:
                        w1 =  ((ye - yd) * (xb - xd) + (xc - xb) * (h - yd)) * r / 2
                        w2 =  (ye - yd) * (xd - xa) * r / 2
                    elif wh > ya and wh <= yd:
                        xf = wh * m
                        xg = xa + (wh - ya) / tan(a2)
                        w1 = ((ye - yd) * (xb - xd) + (xc - xb) * (h - yd)) * r / 2
                        w2 = ((xg - xf) * (wh - ya) * rf + ((xg - xf) * (ye - wh) + (ye - yd) * (xd - xg)) * r) / 2
                    elif wh > yd and wh <= ye:
                        xf = wh * m
                        xg = xd + (wh - yd) / tan(a1)
                        w1 = (((ye - wh) * (xb - xd) + (xc - xb) * (h - wh) + (xg - xd) * (h - wh)) * r + (xg - xd) * (wh - yd) * rf) / 2
                        w2 = ((xd - xf) * (ye - wh) * r + ((xd - xf) * (wh - ya) + (wh - yd) * (xd - xa)) * rf) / 2
                    elif wh > ye:
                        xf = wh * m
                        xg = xd + (wh - yd) / tan(a1)
                        w1 = (((xg - xf) + (xc - xb)) * (h - wh) * r + ((xg - xf) * (wh - ye) + (ye - yd) * (xg - xd)) * rf) / 2
                        w2 = (ye - yd) * (xd - xa) * rf / 2
                    aa1 = w1 * sin(a1) * cos(a1 - a2) + w2 * sin(a2)
                    bb1 = - w1 * cos(a1) * cos(a1 - a2) * 0.74 - w1 * sin(a1) * sin(a1 - a2) * 0.74 - w2 * cos(a2) * 0.74
                    cc1 = w1 * cos(a1) * sin(a1 - a2) * 0.74 * 0.74
                    h1 = (ya + yd + ye) / 3
                    if xc == xb:
                        h2 = (ye + yd + h) / 3
                    else:
                        h2 = (ye + yd + h + h) / 4
                    if h <= 40:
                        alf1 = 1 + h1 / h * (am - 1)
                        alf2 = 1 + h2 / h * (am - 1)
                    else:
                        if h1 <= 0.6 * h:
                            alf1 = 1 + (am - 1) / 3 * h1 / (0.6 * h)
                        else:
                            alf1 = 1 + (am - 1) / 3 + (h1 - 0.6 * h) / (0.4 * h) * (am - 1 - (am - 1) / 3)
                        if h2 <= 0.6 * h:
                            alf2 = 1 + (am - 1) / 3 * h2 / (0.6 * h)
                        else:
                            alf2 = 1 + (am - 1) / 3 + (h2 - 0.6 * h) / (0.4 * h) * (am - 1 - (am - 1) / 3)
                    f1 = ah * 0.25 * w1 * alf1
                    f2 = ah * 0.25 * w2 * alf2
                    aa2 = (w1 * sin(a1) + f1 * cos(a1)) * cos(a1 - a2) + w2 * sin(a2) + f2 * cos(a2)
                    bb2 = - w1 * cos(a1) * cos(a1 - a2) * 0.74 - (w1 * sin(a1) + f1 * cos(a1)) * sin(a1 - a2) * 0.74 - w2 * cos(a2) * 0.74 + f1 * sin(a1) * cos(a1 - a2) * 0.74 + f2 * sin(a2) * 0.74
                    cc2 = w1 * cos(a1) * sin(a1 - a2) * 0.74 * 0.74 - f1 * sin(a1) * sin(a1 - a2) * 0.74 * 0.74
                    dlt1 = bb1 ** 2 - 4 * aa1 * cc1
                    if dlt1 >= 0:
                        k1 = (- bb1 + sqrt(dlt1)) / (2 * aa1)
                        if k1 >= 0 and k1 < kc:
                            kc = k1
                            xd_f = xd
                            yd_f = yd
                    dlt2 = bb2 ** 2 - 4 * aa2 * cc2
                    if dlt2 >= 0:
                        k2 = (- bb2 + sqrt(dlt2)) / (2 * aa2)
                        if k2>=0 and k2 < kc_dz:
                            kc_dz = k2
                            xd_fdz = xd
                            yd_fdz = yd
print('坝高：%6.2f 坝宽：%6.2f 坝坡：%6.2f 水深：%6.2f 湿容重：%6.2f 浮容重：%6.2f 震级：%c' % (h, d, m, wh, r, rf, dz))
print('一般情况下安全系数为 %6.4f 折点坐标 (%6.2f,%6.2f)' % (kc, xd_f, yd_f))
print('地震情况下安全系数为 %6.4f 折点坐标 (%6.2f,%6.2f)' % (kc_dz, xd_fdz, yd_fdz))
f = open('result.txt', 'a')
f.write('坝高：%6.2f 坝宽：%6.2f 坝坡：%6.2f 水深：%6.2f 湿容重：%6.2f 浮容重：%6.2f 震级：%c\n' % (h, d, m, wh, r, rf, dz))
f.write('一般情况下安全系数为 %6.4f 折点坐标 (%6.2f,%6.2f)\n' % (kc, xd_f, yd_f))
f.write('地震情况下安全系数为 %6.4f 折点坐标 (%6.2f,%6.2f)\n' % (kc_dz, xd_fdz, yd_fdz))
f.close
tmp = input('按回车键退出程序')